module mpower
    implicit none
    real(8), private, pointer, dimension(:) :: egv1, egv2 ! , egv1p, egv2p
    integer, private :: size
    integer, private :: fe1=801, fe2=802
    real(8), private :: x1, x2, c1, c2, lk1, lk2
contains

subroutine init_mpower (size_)
    integer, intent(in) :: size_
    size = size_
    allocate(egv1(size))
    allocate(egv2(size))
!    allocate(egv1p(size))
!    allocate(egv2p(size))
!    open(fe1, file='EGV1.txt', status='unknown')
!    open(fe2, file='EGV2.txt', status='unknown')
end subroutine

subroutine add_mpower(egv1n, egv2n, binit)
    real(8), dimension(size) :: egv1n, egv2n
    logical, intent(in) :: binit
    integer :: r1, r2, i
    real(8) :: a, b, c
!    real(8) :: egv1n(size), (size)
    if (binit .eq. .false.) then
        r1 = 0
        r2 = 0
        call find_r1_r2(egv2n,r1,r2)
        a = (egv2n(r1)*egv2(r2) - egv2(r1)*egv2n(r2))
        b = (egv1n(r1)*egv2(r2) + egv2n(r1)*egv1(r2) - ( egv2(r1)*egv1n(r2) + egv1(r1)*egv2n(r2) ) )
        c = (egv1n(r1)*egv1(r2) - egv1(r1)*egv1n(r2))
        x1 = (-b-sqrt(b**2-4*a*c))/2/a
        x2 = (-b+sqrt(b**2-4*a*c))/2/a
        lk1 = (egv1n(r1) + x1*egv2n(r1))/(egv1(r1) + x1*egv2(r1))
        lk2 = (egv1n(r1) + x2*egv2n(r1))/(egv1(r1) + x2*egv2(r1))
        if (lk1<lk2) call swap(x1, x2)
        if (abs(x1)>0.05) x1 = x1/abs(x1)*0.05_8
        egv1 = egv1n + x1*egv2n
        egv2 = egv1n/x2 + egv2n
    else
        egv1 = egv1n
        egv2 = egv2n
    endif
    c1=0;c2=0
    do i=1, size
        c1 = c1 + abs(egv1(i)) !**2
        c2 = c2 + abs(egv2(i)) !**2
    end do
    egv1 = egv1/c1
    egv2 = egv2/c2
!    write(fe1, '(2000(e11.5,1x))') lk1,(egv1(i),i=1,size)
!    write(fe2, '(2000(e11.5,1x))') lk2,(egv2(i),i=1,size)
contains 
subroutine swap(x1, x2)
    real(8), intent(inout) :: x1, x2
    real(8) :: temp
    temp = x1
    x1 = x2
    x2 = temp
end subroutine    
end subroutine

subroutine regen_flux(f1,f2)
    include "xsec.FI"
    include "pspec.FI"
    real(8), pointer, dimension(:,:), intent(inout) :: f1, f2
    real(8) :: f1p(totmesh,ng), f2p(totmesh,ng)
    integer :: i
!    f1p = f1
!    f2p = f2
!    f1 = (f1p + x1*f2p)/c1
!    f2 = (f1p/x2 + f2p)/c2
    
    f1p = (f1 + x1*f2)/c1
    f2p = (f1/x2 + f2)/c2

!    write(fe1, '(2000(e11.5,1x))') lk1,(f1p(i,1),i=1,size)
!    write(fe2, '(2000(e11.5,1x))') lk2,(f2p(i,1),i=1,size)
end subroutine

subroutine find_r1_r2(egv, r1, r2)
    integer,intent(out) :: r1, r2
    real(8),intent(in) :: egv(size)
    integer :: i
    real(8) :: max1, max2
    r1=0
    r2=0
    max1 = 0
    max2 = 0
    do i=1+size*0.1, size*0.9
        if (egv(i)>0) then
            if (egv(i)>max1) then
                max1 = egv(i)
                r1 = i
            endif
        else
            if (-egv(i)>max2) then
                max2 = -egv(i)
                r2 = i
            endif
        endif
    end do
!    do i=1, size/2
!        if (egv(i)>max1) then
!            max1 = egv(i)
!            r1 = i
!        endif
!    end do
!    do i=size/2, size
!        max1 = 0
!        if (egv(i)>max1) then
!            r2 = i
!        endif
!    end do
!    r1 = r1-1
!    r2 = r2-1
end subroutine

subroutine comb(nht, ntrn)
    use define
    use ifport
    use util_mc
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    integer, intent(in) :: nht
    type(particle), dimension(:), intent(inout) :: ntrn
    integer i, l, nq
    real(8) ::xi, x, mu, w1, w2, ws1, ws2
    real(8) :: pi1(size), pi2(size), cj(0:size)
    nq = 0
    do i=1, size
        pi1(i) = abs(egv1(i))
        pi2(i) = abs(egv2(i))
    end do
    pi1 = pi1/sum(pi1)
    pi2 = pi2/sum(pi2)
    cj(0) = 0
    do i=1, size
        cj(i) = cj(i-1) + (pi1(i)+pi2(i))/2
    end do
    
    ws1 = 0
    ws2 = 0
    do i=1, nht
        xi = GetRN()
        do l=1, size-1
            if (xi<cj(l)) then
                exit
            end if
        end do
        x = GetRN() * h(l) + hac(l-1)
        mu= 1- 2*GetRN()
        w1 = pi1(l)*egv1(l)/abs(egv1(l))/(pi1(l)+pi2(l))*2
        w2 = pi2(l)*egv2(l)/abs(egv2(l))/(pi1(l)+pi2(l))*2
        ws1 = ws1 + abs(w1)
        ws2 = ws2 + abs(w2)
        call addntrnm(ntrn, nq, x, l, mu, w1, w2, nn_group())
    end do
end subroutine  

end module
